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We  consider  gravity  waves  at  the  interface  between  two 
uniform,  unbounded  fluids  of  different  densities  in  the 
presence  of  a  current  or  relative  horizontal  velocity  U. 

The  fluids  are  supposed  to  be  immiscible,  incompressible  and 
inviscid,  and  the  motion  is  assumed  to  be  irrotational.  We 
are  concerned  with  the  properties  and  existence  of  finite 
amplitude  two-dimensional,  periodic  waves  of  permanent  form 
which  propagate  steadily  without  change  of  shape.  By  two- 
dimensional,  we  mean  that  the  flow  field  depends  only  on  the 
horizontal  direction  of  propagation,  which  will  be  the  x- 
axis,  and  the  vertical  y-direction.  In  the  field  of  surface 
gravity  waves,  which  is  the  limit  of  the  present  study  when 
the  density  of  the  upper  fluid  is  zero,  it  has  been  found 
recently  that  three-dimensional  waves  of  permanent  form 
exist  and  are  observed  experimentally  (see  e.g.  [5]).  It  is 
expected  that  such  waves  will  also  exist  and  be  important 
for  inter facial  waves,  but  they  will  not  be  considered  in 
the  present  work. 

For  the  purpose  of  calculating  steady  waves,  there  is 
no  loss  of  generality  in  taking  the  speed  of  propagation 
c  parallel  to  the  current  u,  as  an  arbitrary  constant 
transverse  velocity  may  be  linearly  superposed  on  any 
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two-dimensional  steady  wave  without  affecting  its  properties 
(the  stability  characteristics  would,  however,  be 
affected).  The  wave  can  be  reduced  to  rest  by  choosing  a 
frame  of  reference  moving  with  the  wave.  The  problem  is 
then  to  calculate  steady  irrotational  solutions  of  the  Euler 
equations  which  satisfy  continuity  of  pressure  across  a 
common  streamline.  It  follows  from  dimensional  analysis 
that  apart  from  scaling  factors  all  flow  variables  will 
depend  upon  three  dimensionless  parameters: 


h 

L  ' 


pxgL  ' 


(l.D 


where  h  is  the  height  of  the  wave  defined  as  the  vertical 
distance  between  crest  and  trough,  L  is  the  wavelength 
(the  horizontal  distance  over  which  the  flow  field  repeats 
itself  which  in  the  present  work  will  be  the  distance 
between  crests),  p^  and  are  the  densities  of  the 

upper  and  lower  fluid  respectively,  and  g  is  the 
acceleration  due  to  gravity.  For  example,  the  speed  of  the 
waves  is  given  by 


(5L/2.)1/2C(£,  $.) 


(1.2) 


where  C  is  a  dimensionless  function  of  its  arguments.  For 
surface  waves,  where  -  0  and  there  is  dependence  on 
only  one  parameter,  namely  h/L,  it  is  known  that  many 
interesting  and  unexpected  phenomena  exist,  especially  when 
the  wave  steepness  becomes  large.  When  there  is  dependence 
on  three  parameters,  it  is  to  be  expected  that  many  more 
phenomena  are  likely.  However,  in  the  absence  of  exact 
solutions  for  large  h/L,  it  is  a  highly  non-trivial  task 
to  search  a  three-dimensional  parameter  space.  The  results 
to  be  presented  below  are  limited  to  those  phenomena  which 
seem  currently  to  be  of  the  most  interest. 

In  contrast  to  the  voluminous  work  on  surface  waves, 
relatively  little  seems  to  have  been  done  on  interfacial 
waves  of  permanent  form,  and  that  work  seems  to  have  been 
confined  to  the  case  of  zero  current,  i.e.  U  ■  0.  Tsuji 
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and  Nagata  [7]  calculated  Stokes  type  expansions  to  order 
(h/L ) ^,  and  Holyer  [3]  used  the  computer  to  compute  the 
coefficients  in  such  an  expansion  to  order  (h/D37,  and 
then  used  Pad6  approximants  to  estimate  the  behavior  for 
large  h/L.  We  are  not  aware  of  any  work  for  waves  with 
current. 

For  the  mathematical  formulation,  there  is  no  loss  of 
generality  in  taking  g  =  1,  L  =  2it,  and  =  1.  The 
mathematical  problem  is  to  determine  the  x-periodic  velocity 
potentials  and  stream  functions,  for  the  lower  and  upper 
fluid  respectively,  which  satisfy  Laplace's  equation  and  are 
harmonic  conjugate  pairs,  so  that  at  the  unknown  interface 
y  =  Y(x), 

\|>^(x,Y(x) )  =  0,  +2 ( * , Y ( x ) )  *  0  ,  (1.3) 

i-  (V^)2  +  Y(x)  +  b  =  |  p2(V<,2)2  +  pgY ( x )  .  (1.4) 

In  general,  pg  =  p2,  but  we  allow  for  the  possibility  of 
Boussinesq  waves  (in  which  the  inertia  of  the  two  fluids  is 
the  same  and  density  differences  only  matter  when  multiplied 
by  g)  by  setting  p2  =  1  and  Pg  =  0.  Surface  tension  is 
neglected  throughout.  The  quantity  b  is  the  Bernoulli 
constant,  which  by  suitable  choice  of  the  origin  of  pressure 
may  be  set  equal  to  zero  in  the  lower  fluid.  Infinitely  far 
from  the  interface,  we  have 

«*  -Cx,  $2  ~  (U  -  c)x  .  (1.5) 

The  vertical  origin  is  set  by  requiring  that  the  mean 
elevation  of  the  interface  is  zero  and  the  horizontal  origin 
can  be  fixed  by  placing  the  crest  at  x  *  0.  This  problem 
now  appears  to  be  free  of  arbitrary  constants  and  the  wave 
is  determined  by  the  crest  to  trough  height  h.  It  is 
expected  that  isolated  families  of  solutions  exist  in 
connected  regions  in  (h, p2,U)  space,  although  this  does 
not  yet  appear  to  have  been  proved. 

One  question  of  considerable  interest  is  the  domain  of 
parameter  space  in  which  solutions  exist.  Suppose  that  we 
consider  a  fixed  value  of  p2  and  vary  h  and  U.  It  is 
found  that  as  U  increases  with  h  kept  constant  the 
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system  of  equations  describing  steady  solutions  fails  to 
have  a  solution,  even  though  the  'limiting'  wave  profile  is 
smooth  and  exhibits  no  singular  properties.  For  U  >  0 
there  are,  when  solutions  exist,  at  least  two  physically 
distinct  waves  corresponding  to  the  two  wave  speeds  for 
propagation  with  and  against  the  current.  As  U  increases, 
the  wave  propagating  against  the  current  is  'entrained'  by 
the  current  and  at  a  certain  value  of  U,  which  depends 
on  h  and  pj,  the  two  waves  become  identical  and  for 
larger  U  there  are  no  real  solutions  of  the  equations. 
Mathematically,  this  is  like  the  disappearance  of  roots  of  a 
quadratic).  We  shall  term  this  factor  which  limits 
existence  a  'dynamical  limit'. 

The  second  factor  is  what  we  term  a  'geometrical 
limit'.  The  mathematical  formulation  remains  well-behaved 
but  the  solutions  cease  to  make  physical  sense  as  the  wave 
profiles  cross  themselves.  This  occurs  for  fixed  U  and 
increasing  h.  Examples  of  this  phenomenon  are  found  in 
pure  capillary  and  capillary-gravity  waves  [2,1]  for  which 
the  wave  profile  crosses  itself  at  a  critical  value  of  h. 

If  U  *  0,  this  limit  is  going  to  be  different  for  the  two 
solutions  of  waves  moving  with  and  against  the  current.  In 
the  case  of  surface  waves,  this  limit  corresponds  to  a  120* 
cusp.  It  is  easy  to  see  that  except  for  two  special  cases 
(see  §4),  this  cannot  happen  for  interfacial  waves.  Holyer 
[3]  identified  the  geometrical  limit  for  U  ■  0  with  the 
existence  of  a  vertical  tangent.  We  shall  present  evidence 
\  that  waves  can  exist  with  a  vertical  tangent  and  significant 

overhang,  and  the  evidence  indicates  that  the  geometrical 
limit  is  associated  with  the  wave  crossing  itself  when  it  is 
sufficiently  high  for  U  >0. 

2.  WEAKLY  NONLINEAR  WAVES. 

The  properties  of  weakly  nonlinear  steady  waves  may  be 
obtained  by  using  the  Stokes  expansion  in  which  all 
variables  are  expanded  as  power  series  in  h/L.  However, 
the  algebra  can  be  simplified  somewhat  by  using  Whitham's 
variational  approach.  Proceeding  in  the  usual  manner,  one 
finds  after  some  algebra  that  the  average  Lagrangian  is 
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L  "  5  <PB  -  h2  +  a2) 

+  C «*>2  +  P2(Uk-u)2]  -  +  ^g)C“2  +  P2(Uk-w)2] 

-  -jp.  tui2  -  p2(Uk-«)2]  +  0(h6)  (2.1) 

for  the  wave  with  interface  shape 

Y(x)  «*  ^  cos(kx  -  wt)  +  ajCos  2(kx  -  wt)  (2.2) 

The  value  of  a2  is  found  from  3L/ 3a2  =0  to  be 

a2  =  en~-2'p  l[c2  •  p2(u  -  C)21  +  °(h4)  • 

B 

where  C  =  w/k  is  the  phase  speed.  The  dispersion  relation 
for  the  weakly  nonlinear  wave  then  follows  from  3L/ 3h  *  0: 

C2  +  p2(U  -  C)2 

=  (1  -  pb)[1  +  (fP-p-  -  l)2  ♦  jp]  +  0(h4)  .  (2.4) 

For  U  *  0,  the  values  of  C  agree  with  those  in  [7] 

The  values  of  the  energy,  momentum  and  action  densities 
and  fluxes  follow  from  the  expression  (2.1)  for  L  in  the 
usual  way.  In  particular,  the  total  energy  density  E  is 
given  by 

E  *  kCL  -  L  .  (2.5) 

<u 


It  is  to  be  noted  that  for  U  >  0,  the  energy  is  measured 
relative  to  the  energy  of  the  uniform  state  with  a  flat 
interface.  Negative  energies  may  therefore  exist  and  mean 
that  the  energy  of  the  state  with  waves  is  less  than  that  of 
the  undisturbed  flow. 

It  follows  from  the  dispersion  relation  (2.4)  that  for 
linear  waves  (h  ♦  0)  and  given  values  of  p2  and  U, 
there  are  two  solutions  corresponding  to  the  two  roots  of 
the  quadratic  equation  for  C  in  terms  of  p2  and  U.  We 
denote  these  two  solutions  by  C+  and  C_ ,  where 
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C+  >  C_.  For  the  linear  case,  steady  solutions  cease  to 
exist  when  U  exceeds  a  critical  value  UcQ  given  by 

uco  =  Cd  +  P2)U  "  PB)/p2]1/2  (2.6) 


for  which  the  two  wave  speeds  are  equal  with  the  value 

C+  "  C-  “  p2UcO/(1  +  p2}* 

The  values  of  C+  and  C_  are 

P2U  ±  [p2U2  -  (1  +  P2)(P2U2  "  1  +  PbH1/2  •  (2.7) 


For  U  =  0,  the  values  are  equal  and  opposite.  As  U 

increases,  the  speed  of  the  wave  propagating  with  the 

current  originally  increases  but  eventually  decreases. 

The  speed  of  the  wave  propagating  against  the  current 

increases  monotonically  (in  the  algebraic  sense),  becomes 

zero  when  U  =  [(1  -  pb)/p2]  /  and  then  increases  to 

equal  C+  when  U  is  given  by  (2.6).  According  to 

the  linear  approximation,  the  energy  density  E  equals 
12  2 

h  [C  (1+r)  -  rCU],  and  it  is  interesting  that  the  energy 
becomes  negative  when  the  direction  of  the  C_  waves 
changes . 

For  finite  amplitude  waves,  the  two  solutions 
corresponding  to  C+  and  C„  waves  continue  into  two 
families  of  solutions  marked  by  wave  speeds  C+(h, p2,U)  and 
C_(h,p2,U).  For  any  given  value  of  h  and  p2,  there 
will  again  be  a  critical  current  Uc  beyond  which  steady 
solutions  no  longer  exist.  For  the  weakly  nonlinear 
approximation,  this  value  is  given  by 


U 


c 


'W1 


<1  +  p2 } 
4  d  +  p2) 


2 


1/2 

] 


(2.9) 


It  is  noteworthy  that  increasing  h  increases  Uc. 

3.  NUMERICAL  METHODS. 

For  values  of  h  that  are  not  small,  it  is  necessary 
to  employ  numerical  methods.  Three  different  techniques 
were  employed.  The  first  was  to  compute  in  physical  space, 
i.e.  the  interface,  potentials  and  stream  function  were 
expanded  as  Fourier  series  in  x  with  coefficients  which 
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are  exponential  in  y.  The  series  were  truncated  to  N 
modes  and  the  boundary  conditions  were  then  satisfied  at 
N  +  1  equally  horizontally  spaced  points  on  the 
interface.  This  procedure  gives  3N  +  4  equations  for 
3N  +  4  unknowns.  These  equations  were  solved  by  Newton's 
method,  using  continuation  in  either  U  or  h  to  give  the 
first  guesses.  Note  that  this  formulation  is  essentially 
equivalent  to  calculating  numerically  the  coefficients  of 
the  Stokes  expansion  as  done  in  [3]. 

The  second  method  used  the  potential  and  stream 
function  as  the  independent  variables  and  expands  the 
physical  coordinates  as  series  in  these.  The  boundary 
conditions  are  now  satisfied  at  equally  spaced  values  of  the 
velocity  potential  and  the  resulting  system  of  3N  +  3 
equations  in  3N  +  3  variables,  the  expansions  being 
truncated  to  N  modes,  was  also  solved  by  Newton’s  method 
with  continuation  in  U  and  h  employed  to  give  a  first 
guess . 

The  third  method  used  a  vortex  sheet  representation  in 
which  the  unknowns  are  the  shape  of  the  interface  and  the 
dipole  strength  of  the  equivalent  double  layer.  This  gives 
a  nonlinear  integrodif ferential  equation,  which  was  solved 
by  discretization  and  collocation,  the  resulting  system  of 
nonlinear  equations  again  being  solved  by  Newton's  method 
with  continuation. 

For  details,  see  [4,5].  All  methods  worked  extremely 
well  for  small  values  of  h/L,  which  generally  meant 
h  <  0.6,  with  some  dependence  on  Pj  and  U.  (with  our 
scaling,  the  surface  wave  of  greatest  height  has  h  =  0.89). 
The  first  method  was  the  first  to  fail  as  h  increased.  It 
is  of  course  clear  that  this  approach  of  working  in  physical 
space  must  fail  when  the  wave  becomes  very  steep,  but  the 
failure,  marked  by  the  apparent  failure  of  the  Fourier 
series  to  converge,  seemed  to  be  due  to  other  causes.  What 
actually  happened  was  that  the  singularities  of  the  analytic 
continuation  of  the  lower  velocity  potential,  say,  into  the 
upper  half  plane  moved  down  below  the  crest.  In  this  case, 
the  expansion  of  the  velocity  potential  would  have  to 
diverge  near  the  crest,  even  though  the  solution  was 
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perfectly  well  behaved  and  physically  meaningful.  This 
difficulty  would  not  affect  the  other  two  methods  which  were 
used  for  values  of  h  up  to  1.2  for  various  values  of  p 2 
and  U.  For  large  values  of  h,  100  modes  were  used  in  the 
second  method  and  this  seemed  adequate  except  for  the 
largest  h.  The  vortex  sheet  method  with  65  intervals  was 
employed  for  this  case.  This  method  offers  in  principle  the 
advantage  of  being  able  to  concentrate  points  near  regions 
of  high  curvature,  although  this  was  not  done. 

The  accuracy  of  the  calculations  was  checked  by 
comparing  the  results  of  the  somewhat  different  methods  with 
each  other  in  regions  of  apparent  validity  and  by  performing 
the  usual  tests  of  internal  consistency  by  investigating  the 
dependence  on  number  of  retained  modes.  The  calculations 
were  carried  out  on  a  PRIME  750  and  the  CRAY-1  at  NCAR. 

4.  A  SPECIAL  CLASS  OF  SOLUTIONS. 

It  is  interesting  to  note  that  a  special  class  of 
solutions  exist  which  are  simple  transformations  of  the 
well-known  surface  permanent  wave  solutions,  which  have  been 
extensively  studied  both  numerically  and  theoretically  by 
many  authors.  For  each  value  of  p2,  these  solutions 
describe  the  shape  of  the  interface  for  the  C+  case  when 

C+  -  0  -  (1  -  PB)l/2Cs(h)  (4.1) 


where  Ca(h)  is  the  wave  speed  of  the  surface  wave  of 

permanent  frm  for  the  given  wave  height  h.  Since  C+  =  U, 

the  upper  fluid  is  stagnant  in  the  wave-fixed  coordinates. 

The  dynamic  boundary  condition  for  the  motion  in  the  lower 

fluid  then  becomes  that  for  surface  waves  with  a  reduced 

.1/2 


gravity  g(l  -  pb) 


The  velocities  and  wave  speed  are 


therefore  those  of  the  surface  wave  multiplied  by  the  factor 


(1  -  PB) 


1/2 


For  the  C_  branch,  special  solutions  exist  with 
C  -  0,  U  ■  C ( 1  -  pD)/p,]1/2C_(h)  .  (4.2) 


In  this  case,  the  lower  fluid  is  stagnant  and  the  dynamic 
boundary  condition  on  the  motion  of  the  upper  fluid  is  that 
with  a  reduced  upside  down  gravity.  The  wave  profiles  are 
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inverted  surface  waves,  with  negative  gravity  multiplied  by 
the  factor  [(1  -  PgJ/pj]1^2. 

These  special  solutions  have  geometrical  limits  when 
h  =  0.892,  where  the  waves  have  a  corner  at  the  crest  for 
the  C+  wave,  and  a  corner  at  the  trough  for  the  C_  wave, 
with  an  interior  angle  of  120*.  However,  these  special 
geometrical  limits  are  only  for  the  case  when  one  of  the 
fluids  is  moving  with  the  wave.  In  general,  it  is  expected 
(see  below)  that  the  geometrical  limit  is  associated  with 
the  wave  surface  crossing  itself. 

5.  RESULTS . 

The  equations  have  been  solved  numerically  for  various 
values  of  h,  r  =  Pj/p^  an<i  u*  T1le  exi®tence  of  the 
dynamical  limit  was  confirmed,  and  it  was  found  that  the 
weakly  nonlinear  approximation  (2.9)  is  a  good  approximation 
for  values  of  h  up  to  0.6.  A  typical  set  of  results  is 
shown  in  figure  1.  These  results  are  for  r  =  0.5  and 
h  =  0.6,  and  show  the  wave  speeds  C,  total  energies  E, 
and  kinetic  energies  T  for  both  the  +  and  -  waves  as 
functions  of  the  current  velocity  U.  The  existence  of  the 
dynamical  limit  where  C+  *  C_  is  clearly  demonstrated. 

One  feature  of  remarkable  interest  is  the  existence  of 
a  region  of  negative  energy.  This  implies  that  there  will 
be  a  range  of  parameters  in  which  the  energy  of  the  state 
with  finite  amplitude  waves  is  less  than  that  with  the  same 
current  and  a  flat  surface.  Spontaneous  generation  of  such 
flows  is  then  a  definite  possibility.  The  computed  results 
and  linear  analysis  suggest  that  negative  energies  appear 
when  C_  is  zero  and  continue  for  values  of  U  up  to  that 
for  which  the  dynamical  limit  is  reached.  This  aspect  of 
the  solutions  needs  to  be  explored  further  in  detail. 

The  geometrical  limit  or  the  shape  of  the  wave  of 
greatest  height  has  been  addressed  for  the  case  of  zero 
current.  Solutions  were  obtained  for  three  values  of  r 
(1.0,  0.9,  0.1),  using  the  vortex  sheet  method  as  this 
seemed  to  provide  the  best  resolution  when  the  waves  are 
large.  Values  of  the  wave  speed  C  for  +  waves  are  shown 
in  figure  2.  For  the  larger  values  of  r,  it  was  possible 
to  calculate  solutions  with  vertical  slope  and  the  shapes 
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f  steady  interfacial  wave 
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I 


14 


P.  G.  Saffman 


are  shown  in  figures  3  and  4  for  r  =  0.9  and  r  =  1.0, 
respectively.  These  demonstrate  clearly  the  existence  of 
waves  of  permanent  form  with  a  substantial  overhang  region 
in  which  heavy  fluid  lies  on  top  of  light  fluid.  It  is 
interesting  to  note  that  the  fluid  particles  on  the 
interface  between  the  points  of  vertical  tangency  in  the 
overhang  region  are  moving  faster  than  the  wave.  For  the 
smallest  value  of  the  density  ratio,  we  did  not  have 
sufficient  resolution  to  distinguish  the  wave  shape  near  the 
geometrical  limit.  This  difficulty  is  to  be  expected,  since 
the  smaller  the  density  ratio,  the  closer  the  geometrical 
limit  will  be  to  the  120*  cusp  and  the  smaller  the  size  of 
the  overhang  region. 
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